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ABSTRACT 

Using orbital integration and analytical arguments, we have found a new mechanism 
(an "eccentricity trap") to halt type I migration of planets near the inner edge of 
a protoplanetary disk. Because asymmetric eccentricity damping due to disk-planet 
interaction on the innermost planet at the disk edge plays a crucial role in the trap, this 
mechanism requires continuous eccentricity excitation and hence works for a resonantly 
interacting convoy of planets. This trap is so strong that the edge torque exerted 
on the innermost planet can completely halt type I migrations of many outer planets 
through mutual resonant perturbations. Consequently, the convoy stays outside the 
disk edge, as a whole. We have derived semi-analytical formula for the condition for the 
eccentricity trap and predict how many planets are likely to be trapped. We found that 
several planets or more should be trapped by this mechanism in protoplanetary disks 
that have cavities. It can be responsible for the formation of non-resonant, multiple, 
close-in super-Earth systems extending beyond 0.1 AU. Such systems are being revealed 
by radial velocity observations to be quite common around solar-type stars. 



Subject headings: planetary systems: formation - solar system: formation - stars: 
statics 
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INTRODUCTION 



Recent radial velocity surveys indicate that a significant fraction (40-60%) of solar-type stars 
may harbor c lose-in super-Earth s , with masses and p e riods up to ~ 20M m a nd up to a few months, 



respectively (jMavor et al. 



2009 



Bouchv et al 



2009 



Lo curto et al.l |2010| ) . It is suggested that 



in many cases, these super-Earths are members of multiple-planet systems, in which dynamical 
interactions may have influenced their formation although their orbits would not be in mean- 
motion resonances in most systems. Because mm radio observations suggest that the aver aged disk 
mass around solar-type stars is comparable to the mass of the minimum-mass Solar nebula (IHayashi 
198ll ). the total mass of planetesimals would not be enough for in situ accretion of the super-Earths 
in the innermost regions of disks. They may have been accumulated by rocky planets that have 
migrated from outer regions to the proximity of the host stars. Because type I migration is halted 
at the inner disk edge, it may be reasonable that super-Earths accrete near the edge. If the inner 
disk edge is identified as accumulation locations of hot jupiters, its radius may be ~ 0.03-0.05AU. 
However, the observation suggests that averaged orbital radius of the super-Earths is ~ 0.1AU, 
well beyond the disk edge. 



Ogihara &i Idal (J2009J) performed iV-body simulation to study the accretion of planets from 
planetesimals near the disk edge. Inward type I migration of planets is halted either by truncation 
of gas at the disk edge or by resonant perturbation from an inner planet. They found that in 
the case of relatively slow type I migration, 20-30 planets are captured by mutual mean-motion 
resonances and these planets start orbit crossing and giant impacts after disk gas depletion, resulting 
in the formation of several close-in super-Earths. The super-Earths thus formed are kicked out of 
resonances by strong scattering and collisions. This is in contrast to the fast migration case in which 
only several planets survive dur ing the presence of the gas and they remain in stabl e resonant orbits 



even after disk gas is removed (jTerquem Papaloizou 



2007 



Ogihara fc Idall2009h . 



Ogihara &: Idal (|2009l ) also found that in the slow migration case, the innermost planet is pinned 
to the disk inner edge (set at ~ 0.05AU) and the convoy of the resonantly trapped planets extends 
well beyond 0.1AU. Consequently, the super-Earths which form through the giant impacts are 
distributed at ~ 0.05-0.3AU. Thus, the resultant super-Earths are multiple, non-resonant systems 
at ~ 0.1AU, which may be consistent with observed data. 

A question for the result of their TV-body simulation is why the resonantly trapped convoy of 
planets as a whole remain outside the disk edge. Because individual planets (except the innermost 
one) are losing angular momentum through type I migration and the angular momentum is redis- 
tributed throughout the convoy with resonant interactions, large amount of angular momentum 
must be supplied to prevent the planets from pe netrating the disk edge. R everse torque of type I 
migration near the disk edge (IMasset et al.ll2006l ) can halt planets. However. I Ogihara Idal ((2009) 
found the trapping of a convoy outside the disk edge even in the case without the introduction of 
the reverse torque and the effect of the reverse torque does not change the efficiency of the trapping 
and orbital configurations of trapped planets. 
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Here, in order to investigate the trapping mechanism, we simulate in detail the orbital evolution 
of two planets (in some runs, more planets) with various parameters for the disk edge, the planet 
masses, and damping of eccentricities and semimajor axes due to disk-planet interactions. Together 
with analytical arguments, we find that asymmetric damping of eccentricity of the innermost planet 
near the disk edge is responsible for angular momentum input from the disk to the planet with 
such a high rate that it compensates for the angular momentum loss of all the trapped planets due 
to inward type I migration. We call this trapping mechanism an "eccentricity trap." 

Although the eccentricity trap does not occur for some parameter values of the disk and 
planets, once the eccentricity trap occurs, it may be stronger than the trap by the reverse type I 
torque. We derive semi- analytical formula for the conditions necessary to produce an eccentricity 
trap, as a function of planet masses, eccentricity and semimajor axis damping timescales, and disk 
edge width. We also predict how many planets can be trapped. For likely parameter values for the 
innermost regions in protoplanetary disks, the condition for an eccentricity trap may be satisfied, 
while some systems including our solar system might not have had a clear inner cavity in the disk 
and therefore lost the inner planets. 

In section [21 we present analytical argument to clearly show the intrinsic physics of the ec- 
centricity trap. In section [3j we show the results of iV-body simulations. Comparing analytical 
arguments and the numerical results, we derive quantitative conditions for the occurrence of the 
eccentricity trap. In section HI we summarize the results and discuss their implications. 



ECCENTRICITY TRAP 



2.1. Basic Idea 



Consider a planet in an orbit with eccentricity e and semimajor axis a in a gas disk (Fig. [Th). 
We now consid e r changes in a and e ca used by gravitational drag (dynamical friction) from disk gas 
(|Ostrikerlll999l . iTanaka fc Ward! 12004 ; for details, see section [272]) . The tangential velocity of the 
planet near the apocenter in an eccentric orbit is slower than local circular Keplerian velocity when 
the planet moves through an outer disk region. Since the gas velocity is almost equal to circular 
Keplerian velocity, the planet suffers a tailwind. Then, o increases while e decreases, in such a way 
that a(l + e) is almost kept constant (see schematic illustration in Fig. [TJa). On the other hand, 
the planet suffers a headwind when it moves through an inner disk region. Then, both e and a 
decrease such that a(l — e) is almost kept constant (Fig. [lb). Thus, e is decreased by interactions 
with both outer and inner disks, while changes in a due to the the outer and inner disks cancel out. 

If the planet has an eccentric orbit that straddles a relatively sharp disk edge inside which disk 
gas is depleted to create a cavity, it suffers a drag force only when it moves in the outer disk region, 
so that e decreases and a increases on orbit average. Since the specific angular momentum of the 
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planet is 

L = VGM,a(l - e 2 ), (1) 

where M* is the host star mass, the planet apparently gains orbital angular momentum through 
gas drag. 



Note that type I mig ration torque decreases a even in the case of a circular orbit (e.g., IWard 



1986 



Tanaka et al.l 120021 ). In that case, the planet loses L mostly by distant perturbations from 
the outer disk and gains L by those from the inner disk. The loss and gain are almost balanced, 
but the effects of curvature of the system and negative pressure radial gradient slightly enhance 
the loss, resulting in net inward migration (type I migration). Because type I migration is caused 
by residual of opposing torques while the e-damping does not have such cancellation, the timescale 
of type I migration (t a ) is usually much longer than the e-damping timescale (t e ) (see section [3~7Lj) . 

Since the a- increase by the edge effect occurs on a timescale ~ t e , the angular momentum 
gain rate due to the edge effect is much higher than the angular momentum loss rate due to type 
I migration unless e is extremely small (see Eq. [H]). Therefore, if the planet migrates inward 
due to type I migration, it should be "trapped" at the disk edge, as long as orbital eccentricity is 
maintained to some level. If another planet migrates to be trapped in a mean-motion resonance with 
the inner planet trapped at the edge, the resonant perturbations from the outer planet continuously 
pump up eccentricity of the inner one and it enables the eccentricity trap to be maintained. Usually 
other planets also migrate inward and are captured at mutual mean-motion resonances one after 
another to form a convoy of planets. Since the angular momentum gain due to the edge effect is 
significant, it can compensate for the total loss of several planets due to type I migration and trap 
the convoy of planets as a whole outside the edge. 



2.2. Condition for the Trapping 



In this subsection, we derive an approximate condition for the eccentricity trap, in order to 
understand the results of iV-body simulations in section [3l Since in t he iV-body simul a tions, we 
integrate orbits by adding the instantaneous drag forces derived by iTanaka &: Ward! (|2004l ) to 
the equations of motion (rather than secularly changing the orbital elements with orbit-averaged 
torques), we derive the trap condition using their force formulas. 



Goldreich &: Tremaind (|1980i ) , IWardl (j 19881 ) , and lArtvmowica (j 19931 ) derived the damping rate 



of e due to planet-disk interactions. They summed up resonances to obtain an orbit- averaged e- 
damping rate, assuming uniformity of disk gas density on a scale of radial excursion of the planet. 
In the situation we are considering, however, gas density that a planet passes through significantly 
changes during an orbital period, so that we cannot apply their orbit-averaged results. Derivation 
of the formulas for the torque exerted on a planet in an eccentric orbit straddling a disk edge, based 
on full descriptions of Lindblad and corotation resonances, is left for future work. 
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On the other hand, iTanaka fc Ward! (|2004l ) derived instantaneous drag forces caused by grav- 
itational potential perturbations due to density waves at each point of epicycle motion of the 
planet. They showed that the mean e-damping rate (Eq. [5]) obtained by averaging the instanta- 
neous forces over one orbital period ag rees with t hat p redicted by summing up resonant torques in 
the case of locally uniform gas density. lOstrikerl (|1999l ) calculated dynamical friction from gas flow 
onto a particle. Since gas density is higher in downstream than in upstream due to gravitational 
focusing, the particle is pulled back by the gravity o f the higher gas density, which is similar to 
Chand r asekahr ' s dynamical friction in a stellar cluster (jBinnev <fe Tremaindll987l ). ITanaka Ward 
(|2004l ) also showed that spiral waves are always stronger in backside of the epicycle motion (Fig. 1 in 
their paper). The Ostriker's force formula agrees with Tanaka & Ward's one except for a numerical 
factor. 



Because the force formulas by ITanaka &; Ward! (J200J) are local expressions, we adopt them 
both in iV-body simulations and analytic arguments in this paper (see cautions below). In the 
iV-body simulations, the forces are directly included in equations of motion. The formulas for the 
specific damping forces are: 



1 



amp,r 



damp, 



damp, z 



0.78t, 
1 

0.78t, 
1 

0781 



■(2A c r [v e -rn] + A s r v r ) 
{2A c e [v e - rtt] + A s e v r ) 
■(A c z v z + A s z zfl) 



(2) 
(3) 
(4) 



where v r ,ve, and v z are radial, tangential, and vertical components of the planet's velocity, rfi is 
circular velocity of disk gas, which is almost equal to local Keplerian velocity, and the numerical 
factors are given by 





= 0.057 




= 0.176 


A% = 


-0.868 




= 0.325 


K = 


-1.088 


K 


= -0.871 



In the case of lo c ally u niform gas density, integrating -Fdamp.r and i*damp,0 over one orbital period, 
Tanaka h Ward! (|2004l ) derived the averaged e-damping timescale (t e ), 



1 / M 



0.78 VM* 



-i 



-i 



VK/ 



(5) 



where M and M* are masses of the planet and the host star, and S 5 and c s are the gas surface 
density and the sound velocity of the disk gas, respectively. 

These Tanaka & Ward's formulas assumed subsonic flow. We apply the formulas even for 
supersonic cases in which ra dial excursion is wider than the disk edge width that may be comparable 
to disk scale height (H). lOstrikerl (|1999l ) showed that the formula of the drag force must be 
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multiplied by a correction factor, 1/(1 + (f/c s ) 3 ), in the case of supersonic flow. The correction 
factor is consistent with a sup ersonic correction factor for an averaged e-damping rate derived by 



supers 
3), if 



Papaloizou fc Larwoodl (|2000l ). if the relative velocity v is identified by evx- We also carried out 



simulations with the correction factor (section 13, 4p . Tanaka & Ward's formulas also assumed that 
gas density is locally uniform on spatial scale of ~ H. We apply the formulas also for the regions 
near the disk edge in which gas density may vary on a scale of ~ H. The detailed analysis of 
its validity is left for future work, because the purpose of the present paper is to demonstrate a 
potential importance of the eccentricity trap that we firstly found. 

The edge effect described in section 12.11 is quantitatively estimated with local (Hill) approxi- 
mation, as follows. In the Hill coordinates that are comoving with the planet's guiding center at 
a, the tangential velocity of the planet at r = a + x = a(l + ecos(f2ot)) is 

v y = — 2eaQo cos(f2oi) = — 2xOq, (6) 



where f^o = y GM* / a 3 . The Keplerian shear velocity (the local Keplerian circular velocity in the 
Hill coordinates) of gas at x is 

3 

^shear = ~-X^0- (7) 

The tangential relative velocity between the planet and local disk gas is 

Av = Vg - rtl = V y - -Ushear = ~xD,q. (8) 

Since the planet always moves slower than the gas (Av < 0) at x > 0, the drag force accelerates the 
planet's rotation there. At apocenter, the above expression is also easily found in global coordinates. 
The tangential velocity of the planet at the apocenter in the inertial frame is 



L GM,(l-e) 
a(l + e) y a(l + e) 

where L is specific angular momentum of the planet (L = a/ G¥»a(l — e 2 )), while the local circular 
Keplerian velocity is 



Hence, for e <C 1, 



VK ^ = ^Wfe)- (10) 

V0,apo - fK,apo - -^eaO , (11) 



2 

which is identical to Eq. © at the apocenter [x = ea). 

If the disk edge width is sharp enough (ea 3> Ar) and the edge is at x ~ 0, the torque due 
to gravitational drag operates only at x > 0, then the torque due to the e-damping averaged over 
Keplerian time (Tk = 2it/VLq), which we call "edge torque," is 



1 /-Tk/4 

AW k 7jr \ MrF dSimp , e dt, (12) 

J-Tk/4 
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where M is the planet mass. Since vg — rto 
becomes 



(ae/2)r2o cos(riot) and t> r = aef^o sin(f2oi), Eq. 



Then, Eq. {JT2]) is reduced to 



aeto 



o 



0.78*, 



(-.Ag cos(SV) + A| sin(O *)) 



N, 



edge 



M f Tvi/4 A% ean cos(tt t) M ea 2 Q eML 

' adt ~ 



7k 



T K /4 



0.78 



2t P 



where iV e apo is the torque at the apocenter [N t 



e,apo 



2tt tp 27Tt, 



eMa 2 n /2t e ). 



N, 



e,apo 
7T 



(13) 



(14) 



If Ar S> ea, the planet suffers an opposing torque at x < and the net torque must almost 
vanish on orbit average. (In other words, with the integral range from — Tk/2 to Ik/2, the above 
integral vanishes.) Thus, in general cases, the edge torque can be expressed with a reduction factor 
/ (0 < / < 1), which is a function of ea/Ar, as 



N, 



edge 



/ 



eML 
2-kU ' 



(15) 



If the edge is at x > 0, the range of integral in Eq. (JT2J) is —tint to t- mt with t- mt < Ik/4. If the edge 
is at x < 0, ti n t > Ik/4- In either case, the integrated value is smaller than that in Eq. (|14p . which 
means that Eq. (|15p is the maximum value of iV e( j ge . 

In our TV-body simulation, even for planets near th e edge, we apply the conventional inward 
migration with a timescale given by (ITanaka et al.l l2002) 



t a 



a 
a 



1 



2.7+l.lg 



M 



-i 



UK 



to- 



il*) 



where T, g oc r 



Masset et al.l (j20Q6l ) pointed out that near the disk edge, due to the local 



positive surface density gradient, the migration may be reversed to be out ward. The radiative 



effect can also make t he migration ou tward in optically thick disks (e.g., iPaardekooper et al 



2010 



Lyra et al.l |2010| ). Furthermore, iPapaloizou k, Larwoodl (|2000l ) suggested that migration 



is slowed and reversed in supersonic cases (the migration timescale has an additional factor of 
(1 + (evK/c s ) 5 )/(l — (evK I 'cs) 4 )) . Thus, our prescription for type I migration underestimates the 
efficiency of the trapping at the edge. Nevertheless, our iV-body simulation shows that the eccen- 
tricity trap does occur. The effects of the local positive surface density gradient at the disk edge 
and radiation are left to future work. 

In the equilibrium state of the trapping, the edge torque is balanced with the sum of the type 
I migration torque on the planet at the edge and the resonant torque from an outer perturber 
(—N p ). As we assume a conventional type I migration at the edge, the type I migration torque on 



the innermost planet is given by 
at the edge is 



-ML/2t a . Therefore, the condition for the planet to be trapped 
fe ML 



27T t P 



ML 
~2t^ 



N p > 0. 



(17) 
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We first consider the simplest case of two equal-mass (M) planets. In this case, the gravitational 
torque by the inner planet on the outer planet, —N p , balances the type I migration torque on the 
outer planet given by N p ~ ML/2t a . Then, Eq. (fTT|) becomes 

I < £■ (18 > 

Here, we assume that difference in semimajor axes of inner and outer planets is small enough 
and represent them as a mean value, a, in Eq. (|18p for simplicity. Thereby, although we consider 
trapping at the edge, the trapping condition is reduced to that parameterized by the orbit-averaged 
damping timescales of e and a in uniform gas (t e and t a ). 



3. NUMERICAL SIMULATION 

3.1. Comparison between Theory and Simulation 

The above trapping mechanism is quantitatively well demonstrated by simple two-planet sys- 
tems in which the inner planet is located at a sharp disk edge. We have calculated the orbital 
evolution of two planets including effects of type I migration, gravitational drag, and the disk edge. 
The orbits of the planets are integrated by a 4th order Hermite scheme, adding instantaneous drag 
forces Eqs. d2J) , d3j) , and (jH) corresponding to e-damping and a tangential force, 

^mi g ,e = 2^", (19) 



corresponding to a-damping, following lOgihara &: Idal ((2009) 



The parameters in the trapping condition (Eq. |18j ) are t e /t a and Ar/r e dge) where r ec ige is the 
location of the disk edge. If Ar is comparable to the disk scale height (H), 

Ar ' ° s (20) 



^edge V K 

(Note that the disk edge may be truncated by stellar magnetic field and Ar is not necessarily ~ H.) 
From Eqs. ([5]) and (fT6|) . 

t e 2.7 + l.lg f c s \ 2 

— = — < 1- (21) 

t a 0.78 \vk) K 1 

At r ~ 0.05A U corresponding to the inner disk edge, T ~ 1250K in the optically thin limit 



(IHavashil 119811 ). Since v K ^ lSO^/O.OSAU)- 1 / 2 ^^ and c s ~ 2(T/1250K) 1 /2k m /s, c s /v K ^ 0.016 
at r ~ 0.05AU. In an optically thick disk, T can be higher. However, since silicate dust grains 
sublimate and the disk becomes optically thin at T > 1400K, T cannot be higher than that. On 
the other hand, t a can be longer due to non-linear effects. Hence, at r ~ 0.05AU, t e /t a ~ 10 _4 -10~ 3 
and t e ~ IOOTk for M = Af@. We perform calculations with wide ranges of Ar/r e d ge and t e /t a 
with fiducial values of Ar/r e( j ge = 0.01 and t e /t a = 10 -3 . In all runs, t e /^K = 100 and q = 1.5 are 
used, and the gas surface density vanishes with a hyperbolic tangent function of width Ar. 
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The results of the orbital integrations are presented in Figure [21 Two planets with Earth mass 
M ffi are initially placed in almost circular orbits. The top panel shows the orbital evolution of 
the fiducial case with Ar/r e d ge = 0.01 and t e /t a = 10~ 3 . The solid lines represent the semimajor 
axis a, and the dotted lines represent the pericenters a(l — e) and apocenters a(l + e). The inner 
planet migrates to the edge and stays there because T, g vanishes and t a — > 00 there. After that, the 
outer planet migrates to be trapped in the 6:5 mean-motion resonance with the inner one. Once 
the planets are captured in the resonance (at t ~ 45000Tk), the eccentricities of both planets are 
pumped up to e ~ 0.01 and their a's are shifted slightly outward, resulting in an eccentricity trap. 

The eccentricity trap depends on the parameter values. The middle panel shows the result of 
fast migration cases with t e /t a = 10 -2 (Ar/r e d ge = 0.01). The planets are not trapped at the edge 
since the type I migration speed is too fast in this case, although the eccentricities become about 
0.06 at the moment of capture into a resonance. The bottom panel shows the result of smooth edge 
case of Ar/r e d ge = 0.2 (t e jt a = 10~ 3 ). The planets also penetrate into the cavity because the edge 
effect is weak. 

We now derive a trapping condition by using the simulation results. We calculated a torque 
due to the gas drag at each timestep in the orbital integration and found that the orbit-averaged 
torque is well fitted with 

{1 (for er edge > Ar) 

/ eredgc ft a \ ( 22 ) 

y—^T (for er edgc < Ar). 

In the fiducial case with Ar = 0.01r e( j ge , we found that e ~ 0.01, so that / — 1. 

The results with various t e /t a and Ar/r e d ge for two Earth-mass planets are summarized in 
Fig. [3l Crosses represent the cases in which planets are not trapped. Other symbols represent the 
trapped cases. The filled squares, triangles, and circles represent the time-averaged eccentricity of 
the inner planet of e < 0.02, 0.02 < e < 0.03, and e > 0.03, respectively. 

In order to compare the semi-analytical condition, Eq. f)18[) . with the orbital integration re- 
sults, we need a formula for resonantly excited eccentricities. We must also take into account the 
dependence on planetary masses. Let M\ and Mi be the masses of the inner and outer planets, 
respectively. The ec centricity of the inner planet after resonant trapping is expressed in the Hill 



approximation (e.g.. iNakazawa fc Idalll988l ) as 

M 2 _ V2 



- res M 1 + M 2 



(23) 



where e res is the eccentricity of the relative motion (In the Hill approximation, relative motion 
between two Keplerian motions is another Keplerian motion). Since the magnitude of the relative 
eccentricity excited by resonant interaction (e res ) with damping is not theoretically clear, we use 
empirical values. With Mi = Mi = M ffi , we have done the same calculations as Fig. [3] for various 
values of t e [= (10-10 3 )T K ] and t a [= (10 4 -10 6 )T K ]. In Figure HJ e res is plotted against t e and t a . 
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The relative eccentricity can be fitted by 



0-02 ( '^f) V2 • (24) 



The trapping condition for a given mean- motion resonance may be tub < t a , where tub is the 
libration timescale of the eccentricity by the secular perturbation, 

fob (25) 



MJ \ ai J 

where fi = M1M2/ {M\ + M2) is the reduced mass, M and M2 are masses of the inner and outer 



planets, and a± and 02 are their semimajor axes (e.g.. Murray &: Dermottlll999l ). Since tub is larger 
for larger (12/ai, the planets tend to be trapped by a more distant resonance for larger t a . As a 
result, e res is smaller. In section f3.2t we will show through calculations with M\ = M2 = (0.1- 
10)M ffi that the magnitude of e res is almost independent of the planetary mass in the trapped case 
as long as t e /t a is fixed. The fact that all of tub, t e and t a are inversely proportional to planetary 
mass suggests that the resonance at which planets are trapped and hence the values of e res would 
not sensitively depend on planetary masses. However, note that the formula (Eq. [23]) is empirical 
and can strictly be applied only for the range of parameters that we tested. 

Substituting e = e res /2 into Eq. (fT8|) . the trapping condition is reduced to 

for Ar/r cdgc < 0.017z^f, 
'Ar/r edge \~ 2 nni _ 2 (26) 




0.008i/f ( ^l^ ) for Ar/r edge > 0.017^, 



where V2 = 2M2/(Mi + M2) and ^2 = 1 for equal-mass planets. The shaded region in Fig. [3] 
represents the above semi-analytical trapping condition, which is approximately consistent with 
the numerical results of orbital integrations. Note that in iV-body simulations, the magnitude of 
the type I migration torque on the innermost planet should be smaller than ML/2t a (which is 
adopted in Eq. [T7]) because the innermost planet is located near the edge where gas density is 
smaller. Accordingly, the trapping region would be larger than the shaded region. On the other 
hand, neglecting the type I migration torque on the innermost planet in Eq. (117H . we obtain the 
necessary condition for trapping (solid line in Fig. [3j). The transition from the trapping and non- 
trapping should occur in the region between the outer boundary of the shaded region and the 
solid line, which is consistent with the numerical results. Thus, the condition Eq. (I18p is a slightly 
conservative condition. Since the semi-analytical condition includes planetary mass dependence 
through 1/2, we next examine it by comparing orbital integrations and the semi-analytical argument. 



3.2. Dependence on Planetary Mass 



The dependence of trapping condition on the planetary mass is investigated in this subsection. 
Let Mi and M2 be masses of the inner and outer planets. As in the discussion on Eq. (|17p . we 
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represent their semimajor axes (ai and a 2 ) as a mean value, a (~ 01,02), for simplicity. In this 
case, Eq. (fT7|) is 

/d MiL MiL 

where 

M 2 L(a 2 ) _ Mi£(ai) /a 2 y g+1/2 /M 2 \ 2 MjL /M 2 \ 2 

p 2t a (M 2 ,a 2 ) " 2t a (Mi,ai) UJ V^J ~ 2*a(Mi) VAfJ ' 1 J 



With Eq. (|'24p . the trapping condition is given by 

(Mf + MfXMj + M 2 ) 
M\Mi 



< C, (29) 



C = ^ = ^f^V 1/2 . (30) 



where 

, 

TT t e IT V 10 " 3 , 

For a fiducial case with t e /t a = 10~ 3 and Ar/r e d ge = 0.01, e res ~ 0.02 and / ~ 1, so that C ~ 6. 

To confirm the above condition, we carried out additional numerical simulations with various 
planetary masses, (0.1— 10)M®. We used the fiducial parameters, t e /t a = 10~ 3 and Ar/r e d g e = 0.01. 
Since the e-damping time is set such that t e = 100(M/M®) -1 Tk, in order to keep t a /tm, constant, 
we found that the resonance at which planets are trapped is similar and e res ~ 0.02 in most of 
the trapped cases even if M\ and M2 are changed from 0.1M® to lOM®. Note that e res becomes 
larger than 0.02 when Mi ~ CM\ > Mi. Figure [5] summarizes the results on the Mi — M 2 plane. 
The trapped and non-trapped cases are indicated by filled circles and crosses, respectively. The 
shaded region represents the trapping condition (Eq. [29]). The results are consistent with the 
analytical condition except that the trapping is more extended to high M2/-M1 regions, because N- 
body simulations show higher e res than that assumed in the analytical estimate in the high MijM\ 
regions. 

The results might seem rather counter-intuitive. That is, the trap does not occur when the 
outer planet is very small, while the trap does occur even if the outer planet is somewhat more 
massive than the inner one. If the outer planet is too small, the excited magnitude of e\ is too small 
for the edge torque to be effective. On the other hand, even if the outer planet is more massive, 
the trapping can occur, because e% is excited highly enough. For a further massive outer planet, 
the edge torque is no more balanced with angular momentum loss due to type I migration of the 
massive outer planet. Thus, the eccentricity trap works for planets with comparable masses, which 
is often the case during planet formation, because planets actually start migration when migration 
timescale (oc M _1 ) becomes shorter than growth timescale (oc M 1 / 3 ) and hence the migrating 
planets should have similar masses. 
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3.3. Number of Trapped Bodies 

We also investigated systems with more than two planets. The trapping condition becomes 
more severe than in two planet systems, since the edge torque on the innermost planet must balance 
the type I migration torques of all the planets. The number of planets that can be trapped by the 



edge effect increases with increase in the type I migration timescale. In fact. lOgihara Idal ((2009) 



found that 20-30 planetary embryos are lined up by the edge effect in the case of slow migration. 

We estimate the maximum number of trapped planets. For simplicity, we assume equal mass 
bodies. When n planets are trapped near the disk edge, Eq. (f28j) is replaced by 



Tl — 1 

Np 2l7 ML ' (31) 

where we also neglected the a-dependence, although it may not be able to neglected for sufficiently 
large n. Then, Eq. (]17p becomes 



fe TCS 1 _ J_ _ n - 1 
4vr t e 2t a 2t a 

Therefore, 



> 0. (32) 



n< ^TT- (33) 

Through A^-body simulation, we also derived the dependence of n on e res (Eq. [23]) as 



50/ 2 (t e /t a 

n < —JT T7^ ■ ( 35 ) 



Then, Equation (|33|) becomes 

n, < — — I 

10 

To confirm Eq. (|35p , we performed a calculation with n > 3 in the fiducial case with Ar / r e d gc = 
0.01 and t e /t a = 10" 3 . In this case, Eq. ([35]) is reduced to n < 5. Individual planet masses are one 
Earth mass. The orbital evolution is shown in Fig. After the first planet is trapped at the edge, 
subsequent migrating planets are trapped one after another. After 10 5 Tk, four planets are trapped 
by the edge torque. But, when the fifth planet is trapped, the edge torque no longer halts the 
planets outside the edge and the they migrate inward as a whole. The innermost planet is pushed 
into the cavity. After 1.5 x 10 5 Tk, four planets are outside the edge and the edge torque exerted 
on the second planet supports the outer two planets. The result of Af-body simulation is in a good 
agreement with the analytical estimate, n < 5. 

With the value of t e /t a ~ 10~ 5 that Ogihara &; Ida (j2009l ) used in their slow migration case, 



Eq. (|35p shows that n < 500. This is consistent with their result in which they found that all the 
(20-30) migrating planetary embryos are trapped and remain outside the disk edge. 
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3.4. Supersonic Correction 



The Tanaka & Ward's formulas assume subsonic flow. However, for the eccentricity trap to 
be efficient, eccentricity e must be pumped up by resonance to be eo > H, that is, evK > c s . In 



Ostriker 



this c ase, the relative motion between gas and planets can be supersonic. According to 
(j 19991 ). we also carried out runs with the e-damping forces (Eqs. [2]-[4]) multiplied by a supersonic 
correction factor, 1/(1 + (evK/c s ) 3 ). Then, the effective damping timescale i eje ff is defined as 



^e,eff — t e X 



1 + 



C s J 



(36) 



where t e is the Tanaka & Ward's damping timescal e (Eq. [5]). This supersonic correction is consis- 
tent with that obtained by Papaloizou &: Larwood ( 200ol ). 



Because t P . ^ is larger than t e for superson ic cases, this correction is negative for the eccentricity 
trap. However. iPapaloizou &; Larwoodl ([2000) showed that supersonic correction is also applied for 
the a-damping. The correction factor for the a-damping timescale is (l + (evK/c s ) 5 )/(l — (evK/c s ) ), 
which is a stronger function of ev^jcg than that for the e-damping timescale and even changes sign 
to make the migration outward. With the corrections for both a and e damping, the eccentricity 
trap is rather strengthen in supersonic cases. Here, we include the supersonic correction only 
for the e-damping in TV-body simulations, which significantly underestimates the efficiency of the 
eccentricity damping. Nevertheless, the eccentricity trap still occurs for likely values of t e /t a , 
although the region of the eccentricity trap on the Ar/r e dge — t e /t a plane is more restricted. 

Figure [7] summarizes the results with the supersonic correction on the Ar/r e( ^ ge — t e /t a plane. 
(Note that t e which is used in the vertical axis is not the effective damping timescale but that for 
subsonic cases described in Eq. [5].) Figure [7^ and [7b are the cases with c s /vk — 0.03 and 0.02, 
respectively. While the eccentricity trap regions are more restricted, in particular, in Figure 03, 
they still exist. The shaded regions in Fig. [7] represent the (conservative) analytical condition 
(Eq. |18| ) with the effective damping timescale, which are consistent with the numerical results. 



4. CONCLUSION 

We propose a mechanism that we call an "eccentricity trap" of migrating planets and in- 
vestigated the details of the trapping condition both analytically and numerically. This halting 
mechanism is so strong that many resonantly-interacting planets are trapped by the edge torque 
that is exerted only on an innermost planet at the inner disk edge, as long as the edge is sharp 
enough (Ar < r c d gc ; Ar is the edge width and r e d ge is the edge radius.) and the eccentricity 
damping timescale due to planet-disk interaction (t e ) is short enough compared with semimajor 
axis damping timescale (t a ). The explicit trapping condition is given by 

t a 27T 
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where 



/ 



E^edge 

Ar 



(for er edg e > Ar) 
(for er odge < Ar). 



(38) 



The evaluation for the eccentricity of the innermost body, e, is given by Eqs. (J23J) and (J24J). 



As discussed in section [37T] near the inner disk edge (< 0.1AU), the sound velocity scaled by 
Kepler velocity is so small that disk scale height may be sharp (Ar/r e d ge ~ 10~ 2 ) and t e Jt a is 
as small as 1 0~ 4 -10~ 3 even if we assume the formulas for the a and e damping by iTanaka et al 



(|2002l ) and lTanaka fc Wardl (|2004l ). Therefore, the eccentricity trap condition is safely satisfied in 
protoplanetary disks, if they have a cavity. When the motion of the planet is supersonic (evK > c s ), 
damping formulas are modulated. Then, the eccentricity trap is strengthen. 

We also derived a formula for th e maxi mum number of trapped planets outside the edge. It 
explains the result of lOgihara &: Idal (|2009l ) , where 20-30 planets are lined up outside the disk 
edge in the case of slow migration. This is an essential point of lOgihara &: Idal ((2009) 's model 
for formation non-resonant, multiple, close-in super-Earths at ~ 0.1 AU. The innermost planet is 
pinned to the disk inner edge (at ~ 0.05AU) and the convoy of the resonantly trapped planets 
extends well beyond 0.1 AU. After gas dispersal, they undergo orbital crossing and are kicked out of 
the resonances, leading to many giant impacts to increase the planetary mass by more than a factor 
of 10. As a result, a few super-Earths are accreted at ~ 0.05-0.3AU. Since their masses exceed 
several Earth masses after disk dissipation, they avoid the gas accretion necessary to b ecome gas 
giant planets, as we show in a sequential planet formation model in a separate paper (jlda Lin 
2010) . The resultant non-resonant, multiple super-Earths systems at ~ 0.1AU are consistent with 
observed data. The discovered close-in super-Earth systems (GJ 581 and HD 40307) do not have 
any commensurate relationships, and orbits of outer planets extend beyond 0.1 AU, which can be 
formed reasonably by the above model if clear cavities existed in their progenitor disks. 

This trapping mechanism also provides deep insights into satellite form ation around giant 



planets. Based on the slow inflow disk model (e.g., Canup Sz Wardl 2002 . 2006 ). c s /v^ ~ 0.06 and 



we found through a preliminary iV-body simulation that e ~ 0.1 for resonantly trapped satellites. 
Equation ([33]) suggests that a few Galilean satellites are trapped outside the disk edge. If more 
satellites are trapped to violate the trapping condition, the innermost one is pushed into the inner 
cavity. Q Since the disk edge would coincide with corotation radius with the host object's spin, the 
innermost satellite in the cavity loses its orbital angular momentum by tidal torque from Jupiter 
to collide with it. As a result, resonant trapping of a few satellites may be always maintained. 
Because the number of the trapped satellites are relatively small, the satellites may remain stable 
after gas dissipation. An outermost satellite, Callisto, which is not in a resonance, may accrete 



1 In the case of proto-satellite disks, the satellite mass cannot be neglected compared with the disk mass. When 
the total mass of trapped satellites exceeds the disk mass, the trapping may be come no more effe ctive as well. This 
limit for trapping is similar to the above argument of the maximum n ~ a few jSasaki et alJboiol ). 
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i n the outer residual disk. More detailed study of satellite formation is discussed in I Sasaki et al 
()201Cll ) and a forthcoming paper (Ogihara et al. in preparation). 



The validity of applying the e-damping formulas by iTanaka &; Ward! (|2004l ) for nonuniform 
region near the disk edge should b e tested by hydrod ynamical simulations. At the same time, 
reverse torque for type I migration (jMasset et al,ll2006l ). which is not considered here, should also 
be considered as well. Although we simply assumed the disk edge width is similar to the disk 
scale height, disk edge structure may be created by stellar magnetic field. Thus MHD simulations 
are also needed to determine the disk edge structure. Furthermore, the conditions necessary to 
maintain a cavity should also be investigated by MHD simulation. 



ACKNOWLEDGMENT 



We thank fruitful and stimulating discussions with Steve Lubow, Bill Ward, Doug Lin, John 
Papaloizou, Aurelien Crida, Clement Baruteau, and Sijme-Jan Paardekooper during our participa- 
tion in "Dynamics of Discs and Planets" workshop at the Newton Institute at Cambridge University, 
where some of this work was done. We also thank the anonymous referee for helpful comments, and 
thank Hidekazu Tanaka, Taku Takeuchi and Takayuki Muto for valuable suggestions. This work 
was supported by Grant-in- Aid for JSPS Fellows (2008528). 



-16- 



REFERENCES 

Artymowicz, P. 1993, ApJ, 419, 166 

Binney, J., k Tremaine, S. 1987, Nature, 326, 219 

Bouchy, F., Mayor, M., Lovis, C, Udry, S., Benz, W., Bertaux, J. -L., Delfosse, X., Mordasini, C, 
Pepe, F., Queloz, D., k Segransan, D. 2009, A&A, 496, 527 

Canup, R. M., & Ward, W. R. 2002, ApJ, 124, 3404 

Canup, R. M., k Ward, W. R. 2006, Nature, 441, 834 

Goldreich, P., k Tremaine, S. 1980, ApJ, 241, 425 

Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 

Ida, S. k Lin, D. N. C. 2008, ApJ, 673, 484 

Ida, S. k Lin, D. N. C. 2010, ApJ, 719, 810 

Lo Curto, G., Mayor, M., Benz, W., Bouchy, F., Lovis, C, Moutou, C, Naef, D., Pepe, F., Queloz, 
D., Santos, N. C, Segransan, D., k Udry, S. 2010, A&A, 512, A48 

Lyra, W., Paardekooper, S. -J., k Mac Low, M. -M. 2010 ApJ, 715, L68 

Masset, F. S., D'Angelo, G., k Kley, W. 2006, ApJ, 652, 730 

Mayor, M., Udry, S., Lovis, C., Pepe, F., Queloz, D., Benz, W., Bertaux, J. -L., Bouchy, F., 
Mordasini, C, k Segransan, D. 2009, A&A, 493, 636 

Murray, C. D., k Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge Univ. 
Press) 

Nakazawa, K., k Ida, S. 1988, Prog. Theor. Phys. Suppl., 96, 167 
Ogihara, M., k Ida, S. 2009, ApJ, 699, 824 
Ostriker, E. C. 1999, ApJ, 513, 252 

Paardekooper, S. -J., Baruteau, C, Crida, A., k Kley, W. 2010 MNRAS, 401, 1950 
Papaloizou, J. C. B., k Larwood, J. D. 2000 MNRAS, 315, 823 
Sasaki, T., Ida, S., k Stewart, G. R. 2010, ApJ, 714, 1052 
Tanaka, H., Takeuchi, T., k Ward, W. R. 2002, ApJ, 565, 1257 
Tanaka, H., k Ward, W. R. 2004, ApJ, 602, 388 



-17- 



Terquem, C, & Papaloizou, J. C. B. 2007, ApJ, 654, 1110 
Ward, W. R. 1986, Icarus, 67, 164 
Ward, W. R. 1988, Icarus, 73, 330 



This preprint was prepared with the AAS IATjrjX macros v5.2. 



-18- 




Fig. 1. — Schematic illustration of expansion/shrink of a due to e-damping. (a) The solid line 
represents an eccentric orbit with semimajor axis a around a star. The shaded toroidal region 
represents a gas disk. The guiding center, which divides the disk into two areas, is expressed by the 
dashed line. At apocenter, the tangential velocity of the planet is slower than the local gas velocity, 
while the orbital velocity is faster than the gas velocity at pericenter. (b) The planet suffers a 
tailwind when it is moving in the outer disk, which increases a while decreases e with the apocenter 
a(l + e) kept almost constant, (c) In the inner disk region, the planet suffers a headwind, leading 
to a decrease of a, with the pericenter o(l — e) kept almost constant. 



J I L 



2xl0 4 4xl0 4 6xl0 4 8xl0 4 10 5 




2xl0 4 4xl0 4 6xl0 4 8xl0 4 10 5 

tlT* 

Fig. 2. — Orbital evolution of the planets. Solid lines are semimajor axes, and dotted lines are 
apocenters and pericenters. Top: The result of the fiducial case with Ar/r e d ge = 0.01 and t e /t a = 
1CT 3 . Middle: The result of a fast migration case with t e /t a = 10 2 . Bottom: The result of a 
smooth edge case with Ar/r e d ge = 0.2. 
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Fig. 3. — The results of trapping for various values of Ar/r e d ge and t e /t a . Crosses indicate that 
planets are not trapped, while other filled symbols indicate that planets are trapped. The filled 
squares, triangles, and circles represent the time-averaged eccentricity e < 0.02, 0.02 < e < 0.03, 
and e > 0.03, respectively. The shaded region represents the theoretically predicted trapping regions 
(Eq. [26]). The solid line represents the upper limit (necessary condition) for the eccentricity trap. 
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Fig. 4. — The relative eccentricity excited by resonant trapping is plotted against t e and t a through 
numerical simulations. In these calculations, M\ = M2 = M® is assumed. The contours are drawn 
from e res = 10 -2 5 to 10 -1 ' 2 with intervals Alog 10 e res = 0.1. 




Fig. 5. — The results of trapping for various masses of inner and outer planets (Mi, M2). Filled 
circles represent the cases where planets are trapped, while crosses represent non-trapped cases. 
The shaded region represents the theoretically predicted trapping regions given by Eq. (I29p . 
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Fig. 6. — Orbital evolution of multiple (n > 2) Earth-mass planets in the fiducial case (Ar/r e d ge = 
0.01 and t e /t a = 10~ 3 ). Solid lines are semimajor axes, and dotted lines are pericenters and 
apocenters. 
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Fig. 7. — Same as Fig. [3l but with supersonic correction factor, (a) Case with c s /vk — 0.03. (b) 
Case with c s /vk — 0.02. The shaded region is the analytically predicted trapping regions. Note 
that t e is not the effective damping timescale. 



